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Abstract 

Molecular Dynamics and Statistical Mechanics make possible a particle-based understanding of 
Thermodynamics and Hydrodynamics, including the fascinating Loschmidt contradiction between 
time-reversible atomistic mechanics and the time-irreversible thermodynamic dissipation incorpo- 
rated into macroscopic fluid and solid mechanics. 
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I. INTRODUCTION 



Among the many problem areas in mechanics the study of instabilities and irreversible 
processes seem particularly interesting. Engineering mechanics exists as a discipline because 
failures of structures cost so many lives. The analysis of local Lyapunov instability gives a 
means for localizing and predicting catastrophic failures so that there is a decidedly practical 
engineering aspect of this fascinating scientific research area. 

The goal we pursue here is to develop models that help us to understand. It is pro- 
foundly interesting that the small scale microscopic models of material behavior (ordinary 
Newtonian mechanics) are time-reversible while the macroscopic models of the same thing 
(finite-element and finite-difference fluid mechanics and solid mechanics) are irreversible. 
The tools from nonlinear dynamics and chaos are useful in analyzing these two kinds of 
description. 

Here we describe the basic building blocks for particle simulation and point out the ways 
that these time-reversible simulations already lead to time-irreversible behavior. Most of 
the examples treated here are also described in our three books on computational statistical 
mechanics, smooth-particle applied mechanics, and time reversibility, computer simulation, 
algorithms, and chaosi - -. 

II. ALGORITHM FOR CONSERVATIVE PARTICLE MECHANICS 

No special tricks are necessary to get started with particle-based simulations. Microscopic 
mechanics can provide us with accurate particle trajectories { q(t) } . All we need to do is 
to integrate Newton's ordinary differential equations of motion , 



Here $ is the potential energy, a function of the coordinates { q } . Alternatively, we 
can obtain an equivalent coordinate- momentum description { q(t),p(t) } by integrating 
Hamilton's first-order ordinary differential equations : 



Both these approaches are time-reversible. That is, a movie of the motion, played backwards, 
satisfies exactly the same equations (with the values of q and p in reversed order and with 



{ mq = F(q) 



V$ } . 



{ q = +(&H/dp) ; p 



(dH/dq) } . 
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the sign of p changed also). A movie is an excellent analog of numerical simulation. Both 
the simulation and the movie are sets of discrete records of coordinates at discrete values of 
the time, separated by the "timestep" At . In addition to the basic algorithm keep in mind 
that three crucial questions remain to be answered: [1] what are the initial conditions, [2] 
what are the boundary conditions, and (most important of all) [3] what is the problem to 
be solved? 

Macroscopic continuum mechanics is based on the three conservation laws for mass, 
momentum, and energy : 

p = - p V • u ; pit = -V • P ; pe = -Vit : P - V • Q . 

Here p is density, u is velocity, e is energy per unit mass, P is the pressure tensor (force per 
unit area, necessarily a symmetric second-rank tensor) , and Q is the heat flux vector (energy 
flow per unit area). All these variables are continuous functions of space and time. Both 
P and Q , as well as p , u , and e , are defined in the comoving frame, a coordinate frame 
moving with the local velocity u(r,t) . Finite-difference approximations to the gradients on 
the righthand sides of the three continuum equations, evaluated at a discrete set of spatial 
mesh points, reduce the partial differential equations to ordinary ones, which can then be 
solved with Runge-Kutta integration. Again, the hard part of the problem is the same: 
what to do and how to implement the initial and boundary conditions. 

Different materials can be described by different types of constitutive relations (elastic, 
plastic, viscous, ...) giving P and Q in terms of the basic { p,u,e } set, together with 
their time derivatives and spatial gradients. Time-reversed movies of solved macroscopic 
problems look "funny" and make no sense. This is because the underlying phenomenolog- 
ical constitutive relations are typically irreversible. The simplest most familiar irreversible 
examples are Newtonian viscosity and Fourier heat conduction : 

p = [ p cq - AV • u ]I - r)[ Vu + W ] ; Q = -kVT . 

In the symmetrized velocity gradient Vw* is the transpose of Vw . / is the unit tensor. 
A boxed conducting fluid, with that fluid initially in motion, (a Rayleigh-Benard flow, 
for instance, but with the box suddenly insulated and with the accelerating gravitational 
field suddenly switched off) described with a shear viscosity r\ and a heat conductivity k 
eventually comes to an isothermal state of rest. Evidently the reversed movie of this decay 
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makes no sense and would correspond to an illegal "something from nothing" contradicting 
the Second Law of Thermodynamics. Of course, with the right initial conditions and the 
right boundary conditions, one can indeed observe tornadoes! For an Lx L (two-dimensional) 
system, with kinematic viscosity and thermal diffusivity of order D the initial gradients decay 
exponentially, ~ e -Dt/L 2 _ rj^g reversed movie, with its exponential growth, ~ e +Dt/L 2 ^ j g 
simply wrong. 

In applications of mechanics to simulation we strongly recommend the use of the fourth- 
order Runge-Kutta integrator because it is easy to use and to modify for the treatment of 
open systems interacting with their environments. If the focus is on the time-reversibility of 
conservative Newtonian systems it is useful to consider a very simple, yet rigorously time- 
reversible, integrator discovered by Levesque and VerleiA In order better to understand 
the coexistence of the reversible microscopic and irreversible macroscopic views we adopt 
Levesque and Verlet's "bit-reversible" algorithm. This approach generates a numerical tra- 
jectory in an integer coordinate space, by rounding off the acceleration terms : 

{ q n+1 - 2q n + q n _i = [ (F (q n ) (At) 2 / 'm ] in te ge r } • 

The subscripts indicate the time, in units of the (integer) timestep At . The initial conditions 
to start this algorithm are the coordinates at two successive times. 

A simple illustration of the algorithm follows a harmonic oscillator trajectory, using unit 
mass, force constant, and timestep At : 

q+ - 2g + q- = -qo — ► q+ = go - q- ■ 

The solution of repeating coordinates { +1, +1, 0, — 1,-1,0,...} , is typical, and illustrates 
the fact that no matter what the initial conditions, the solution is both periodic (for chaotic 
problems, the length of the period is of order the square root of the number of states) 
and reversible. The algorithm is a faithful analog of classical deterministic time-reversible 
mechanics. If momenta are desired they too can be approximated accurately from the 
coordinate values : 



Po 



(q+ - Q-) 



(q++ - q—) 

AAt 



2At 

Figure 1 compares the energy calculated with these momenta with calculations based on 
the Beeman and velo city- Ver let algorithms. Our formulation (small dots in the Figure) is 
clearly an improvement. It is evident that a promising research area lies in the development 
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FIG. 1: The Beeman and Velocity-Verlet approximations to the harmonic oscillator momentum are 
compared with the more nearly accurate formula given in the text (at the left). The corresponding 
total energies are shown at the right over two oscillations with At = 1 . The approximation 
mentioned in the text is the best of the three and corresponds to the smallest dots. 



of higher-order bit-reversible algorithms combining coordinates, velocities, and accelerations 
from more than three successive times. 



III. IRREVERSIBILITY FROM TIME-REVERSIBLE MOTION EQUATIONS? 

Is there any chance of detecting irreversibility with such a time-reversible algorithm? 
Oddly enough, there is! It is based on the analysis of Lyapunov instability, looking in the 
neighborhood of the trajectory, not just at the trajectory itself. Such a nonlocal analysis 
necessarily depends upon the imagination and the information contributed by an observer 
of the motion. Let us take an example, a maximally irreversible situation described by 
time-reversible, even bit-reversible, Newtonian mechanics. Consider the pair of Shockwaves 
launched by the collision of two mirror-image fluid samples. See Figure 2 for four snapshots 
of such a problem. Initially the velocities are ±u . Eventually the periodic L x L system 
shows no more systematic motion - the initial kinetic energy has been completely converted 
to internal energy (heat) : 

(v?/2) -> e . 

Just as before, any portion of the developing trajectory can be reversed precisely and exactly 
despite the Lyapunov instability (exponential growth of perturbations) of the dynamics. 
There is a vast literature^- on the quantification of Lyapunov instability, the exponen- 
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FIG. 2: Snapshots from the twofold compression of a 40 x 40 cold crystal with unit density and 
with the pair potential (ft = (1 — r 2 ) 4 . Initially half the particles move to the right and half to the 
left, at speed 0.875 . The second snapshot, time = 11.5, and located at bottom right, corresponds 
to the maximum twofold compression. Such "irreversible" motions can be reversed precisely, for 
as long as desired, with the Levesque-Verlet bit-reversible integration algorithm. 

tially sensitive deformation of comoving hypervolumes in q space, p space, or { q,p, } phase 
space. For N particles in two dimensions the 4iV-dimensional phase-space motion defines 4iV 
local Lyapunov exponents. The fact that these "local" exponents depend upon the chosen 
coordinate system can be viewed as a disadvantage or as an opportunity. Again, there are 
many promising research problems suggested by this observation. Optimizing the analysis 
is certainly a useful and stimulating activity. 

The largest Lyapunov exponent - we will call it X\(t) - associated with the motion can be 
found by following two nearby trajectories in time. The primary or "reference" trajectory 
can be generated with bit-reversible dynamics, so that it is possible to extend it as far 
as desired into the future or the past^. The dynamics of a nearby "satellite" trajectory 
is restricted by constraining the satellite trajectory to stay within a fixed distance of the 
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reference trajectory. The satellite dynamics can readily be generated with Runge-Kutta 
integration, rescaling the separation between the two trajectories at the end of each time 
step. The local Lyapunov exponent is : 

(Ar) before 



Ai(t) = 



At 



In 



(At] 



after 



, Ar — ('"satellite '"reference! 



Although the motion equations are perfectly reversible for both the reference and the satel- 
lite, the reversed satellite trajectory turns out to be totally unlike the forward one, if the 
system is a nonequilibrium system. Both the local Lyapunov exponent associated with the 
instability and the identities of those particles making above-average contributions to the 
offset vector separating the trajectories , Ar = r sate iiite — '"reference, are qualitatively different. 
In a typical shockwave simulation of the type shown in Figure 2 the number of these more- 
heavily-weighted particles is about twice as great in the reversed motion as in the forward 
one. 

This is an extremely interesting result. No doubt it suggests various "Arrows of Time" 
which can be constructed based on the structure of nearby trajectories (which react to the 
past, not the future). A study of irreversible flows from this standpoint should shed light 
on the reversibility paradox for simple Newtonian and Hamiltonian systems. 



IV. IRREVERSIBILITY FOR TIME-REVERSIBLE "OPEN" SYSTEMS 

In order to control nonequilibrium states, in particular nonequilibrium steady states, 
it is necessary to do work and/or to exchange heat, with the system of interest. The 
dynamics becomes a little more complicated due to these interactions, but the interpretation 
compensates by becoming simpler. Here we take up the description of "open systems" . 



A. Microscopic Pressure, Heat Flux, and Temperature 

"Open" systems have mechanical and/or thermal connections to their environment, open- 
ing up the possibility of simulating processes including thermodynamic work and the flow 
of heat. Analysis of these systems requires microscopic analogs for all the continuum vari- 
ables. Density, velocity, and energy are the simplest of these. We adopt the smooth-particle 
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averaging method discovered by Lucy and Monaghan in 1977—^ : 

P( r >^) = E m * w ( r_r *) ; p(r,t)u(r,t) = ^mj^w(r-rj) ; p(r,t)e(r,t) = ^rriieiw(r-ri) . 

i i i 

The particle energies { e« } are necessarily defined in the comoving frame, which moves with 
velocity u(r, t) . A useful weight function, with a range h which can be optimized, is Lucy's, 
here normalized for two space dimensions : 

w (r) = (5/7r/i 2 )(l + 3z)(l - zf z = \r\/h . 

This weight-function approach guarantees the continuity of the first and second spatial 
derivatives of the field variables and lends itself to optimization studies. 

In addition to the basic mass, momentum, and energy, several other variables need to 
be considered in order to compare microscopic and macroscopic simulations. Unlike their 
continuum cousins, the pressure tensor and the heat flux vector from molecular dynamics 
are respectively even and odd functions of the time : 

pv = E (rF)o- + E(pp/™)< ; Qv = E ml F n • (p< + ^)/ 2 ] + E( e ?/™)* • 

pairs i pairs i 

Here = r, — rj and Fij is the force on Particle i due to its interaction with Particle j. 
The individual particle energies { } include half of each particle's pair interactions with 
its neighbors. 

Temperature needs a definition too. The usual equilibrium definition, based on entropy, is 
useless away from equilibrium where entropy has no consistent definitioni>2. At equilibrium 
Temperature can be defined in many ways, all based on Gibbs' statistical mechanics or 
Maxwell and Boltzmann's kinetic theory. The even moments of the velocity distribution 
are examples. In addition to these there are also configurational definitions. The simplest 
"configurational temperature" is based on an identity from Landau and Lifshitz' text^ 1 ^ : 

This definition follows from an integration by parts in Gibbs' canonical ensemble. If the 
differentiation indicated by V is carried out in momentum space the Landau-Lifshitz formula 
gives the usual kinetic-theory definition of temperature , mkT xx = ( p x ) . If instead the 
gradient is carried out in coordinate space the "configurational temperature" depends on 
the first and second derivatives of the potential function governing the motion : 

^^configurational = ( F 2 )/( V 2 $ ) . 
8 



One-body or many-body configurational temperatures, either scalar or tensor, can be defined 
in this way. But an evaluation of them for the shockwave problem reveals divergences. 
Typically particle values of V 2 $ frequently alternate between positive or negative values, 
so that the corresponding configurational temperatures frequently diverge! Configurational 
temperature also has unphysical undesirable contributions arising from rotation whenever 
Coriolis' or centrifugal forces are significant. 

The simplest definition for temperature is the kinetic second-moment one. It is based on 
a mechanical model of a working ideal-gas thermometer. In that instance a relatively heavy 
mass-M "system atom" interacts with a collection of light-weight mass-m "ideal-gas ther- 
mometer" particles characterized by an unchanging equilibrium Maxwell-Boltzmann distri- 
bution with temperature T . Kinetic theory shows that the averaged effect of such collisions 
causes the system-atom velocity to decay while its mean-squared velocity approaches the 
equilibrium value for the temperature T : 

( v x } oc ~{v x /t) ; ( v x v x ) oc [ (kT xx /m) - v 2 x )/r ; [ for m << M ] . 

Accordingly, we adopt the kinetic definition of temperature in what follows. With temper- 
ature defined we can proceed to devise "thermostats" able to control it. 

B. Time-Reversible Deterministic Thermostats 

The first of the deterministic mechanical thermostats was Woodcock's isokinetic 
thermostat^, implemented by rescaling the velocities at the end of each timestep. Much 
later it was discoverecr^ 1 ^ that a continuous time-reversible version of this thermostat could 
be implemented with a time-reversible friction coefficient ( : 

v = F( q ) - Cp ; C = EOF ■ v )l J> 2 /™) (d/dt) £(p 2 /™) = o . 

This "isokinetic" thermostat can be applied to one or many particles and to one or many 
space directions. 

An illustrative application is the "Galton Board"- 1 ^ 1 ^, in which a single particle is accel- 
erated through a lattice of scatterers but constrained to move at constant speed. Overall, 
the potential energy drops. Because the mean value of the friction coefficient is necessarily 
positive, the phase-space probability density collapses onto a multifractal strange attractor, 
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quantifying the rarity of nonequilibrium phase-space states. This approach to tempera- 
ture control is often termed the "Gaussian" thermostat because Gauss' Principle (of Least 
Constraint) gives this thermostat when applied to the problem of constraining the kinetic 
energy^. Reference 15 is a detailed discussion of the model (summarized in References 1 
and 2). This work clearly shows the fractal nature of the phase space (with vanishing phase 
volume) that results when the dynamics is thermostated. Fancier thermostats, based on 
statistical mechanics, can be found. 

In 1984 Shuichi Nose discovered a precursor of the best of them, a thermostat^ with a 
more elaborate basis in Lagrangian and Hamiltonian mechanics, but somewhat disfigured by 
a novel "time-scaling variable" s = dt \d/dt new . His thermostat imposed Gibbs' canonical 
phase space distribution at equilibrium rather than the less-usual isokinetic one. A sim- 
plification of his equation of motion, without the useless time-scaling, likewise contained a 
friction coefficient, which itself obeyed an evolution equation depending upon past values of 
the kinetic energy : 

{ V = F(q) - (p } ; C = [ (T({p})/T ) - 1 }/r 2 [ Nose - Hoover ] . 

The relaxation time r is a free parameter determining the time required for the thermostat 
forces {— (p} to bring the kinetic temperature T({p}) to the desired thermostat temperature 
T . Just as in the isokinetic case the nonequilibrium averaged friction coefficient for this 
"Nose-Hoover" mechanics is positive, leading once again to multifractal phase space distri- 
butions away from equilibrium. There is an extensive somewhat mathematical literature 
having to do with picking the "right" relaxation time or the "right" thermostat. 

The Gaussian and Nose-Hoover thermostats are particularly useful for controlling 
nonequilibrium problems, such as shear flows and heat flows, and the Rayleigh-Benard 
problem combining them. The definition of temperature depends upon the definition of 
the local velocity u(r,t) . A straightforward definition of velocity, which nicely satisfies the 
continuity equation exactly^, can be based on smooth-particle weighting functions : 

u(r,t) = ^Viw{\r - Ti\) /Y^,w{\r - Ti\) . 

i i 

The Gauss, Nose, and Nose-Hoover thermostats can all be related to Hamiltonian me- 
chanics. Dettmann, together with Morriss, carried out much of this work, with later con- 
tributions by Bond, Laird, and Leimkuhler, and then by Campisi, Hanggi, Talkner, and 
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Zhai>i£~— . All of them helped to clarify the connections of time-reversible thermostats with 
standard Hamiltonian mechanics. This work leads to the conclusion that many different 
thermostats can be used at equilibrium but that some of them fail in nonequilibrium situ- 
ations, even in situations close to equilibrium. Just as in real life the failures, rather than 
the successes, are the more newsworthy subjects. Let us turn to some examples. 

C. Thermostat Failures — Oscillators, Heat Conduction, and the cf) 4 Model 

A very stimulating "log-thermostat" has just been described by Campisi, Hanggi, Talkner, 
and Zhan^i They pointed out that the microcanonical (constant energy) ensemble distribu- 
tion for a logarithmic potential generates (at least formally) the Maxwell-Boltzmann velocity 
distribution : 

r+oo 

= kT In q -> / dq5\ 2H - kT In q 2 - (p 2 /m) ] oc exp[ (H Q /kT) - {p 2 /2mkT) ] 
Jo 

Because the dynamics of this thermostat is unstable, there being nothing to keep q away 
from the origin, in applications they recommend using kT ln(q 2 + 5 2 ) , where 5 is sufficiently 
small. 

Our effort to use this thermostat for a nonequilibrium heat flow problem failed. Connect- 
ing a cold and a hot log-thermostat to opposite ends of a two-particle 4 chain gave different 
temperatures at the two ends, but no heat flux at all. The problem is that the Hamil- 
tonian log-thermostat is unable to replicate the phase-space contraction associated with 
dissipative systems. There are some other examples of such failures. Leete and Hoover's 
HamiltoniarJ2i22 , 

Uul = sl±K{p)K + <%) - K , 

keeps the kinetic energy, ^{mq 2 /2) constant, equal to K Q . The configurational temperature 
can alternatively be kept constant using a special Hamiltonian. In both these cases a cold 
and a hot thermostated region, in contact with Newtonian regions, gives no heat flux at all 
despite huge temperature differences. The lesson is that Hamiltonian mechanics is not able 
to describe dissipation properly. 
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V. THERMOSTAT SUCCESSES" OSCILLATORS AND COMPLEX SYSTEMS 



A "good" thermostat should, for instance, be able to provide good solutions of the 
Rayleigh-Benard problem, heat transfer through a compressible fluid in a gravitational field. 
It should also be useful in treating small-scale "toy problems" . The simplest thermostat test 
problems are [1] a harmonic oscillator— with a coordinate-dependent temperature : 

T(q) = 1 + etanh(g) ; 

and [2] the flow of heat through a 4 chain of particles : 

« = (P?/2m) + W/4) ] + J>y , 

i i<j 

where is a nearest-neighbor Hooke's-Law potential and where the first few and last few 
particles in the chain are thermostated with a Gaussian or Nose-Hoover or another thermo- 
stat. This model is a specially good one to study because it is known to satisfy Fourier's 
law, even in one dimension. A comparison of seven different thermostat methods showed 
that the <j) 4 problem is well-posed and relatively easy to solve 25 . 

Some Hamiltonian-based thermostats are ineffective for nonequilibrium problems 2 ^ - — and 
it is useful to understand why. At equilibrium a given temperature and volume imply corre- 
sponding values of the kinetic and potential energies. This is also true for particular states 
away from equilibrium, even where there is no longer a unique equation-of-state relation. 
Using a Hamiltonian thermostat away from equilibrium one can independently specify the 
kinetic energy and the potential energy or the temperature and the energy. This additional 
freedom contradicts the notion of thermodynamic state and can lead to very strange results^. 
Constraining the configurational temperature or using a version of Hamiltonian mechanics 
to constrain the kinetic energy discovered by Hoover and Leete provide temperature pro- 
files that make no sense. The log-thermostat is another demonstration that Hamiltonian 
mechanics is a poor choice for thermostats. This paradoxical situation is the symptom 
of two incompatible requirements on the dynamics: [1] Liouville's Theorem requires that 
the phase-space motion be incompressible; [2] Heat flow consistent with the Second Law of 
Thermodynamics requires that the phase volume decrease to zero. 

Several of the thermostats have no problem with generating heat flows and solve the 
problem of decreasing phase volume by generating strange attractors in the phase space. 
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0.4 < e < 1.0 




-4 < p < +2 

FIG. 3: £(p) is plotted for fifteen equally-spaced values (e = 0.44 to 1.00) of the maximum tem- 
perature gradient for the (q,p, C) nonequilibrium Nose-Hoover oscillator. All these data correspond 
to fully-converged limit cycles. 

Let us consider what might appear to be the simplest of these problems, the Nose-Hoover 
oscillator 24 with a temperature gradient 1 - : 

q — p ; p — — q — Qp ; £ = p 2 — T(q) ; T(q) = 1 + etanh(g) ; < e < 1 . 

For e > 0.4 the motion is a one- dimensional limit cycle with ( ( ) positive. The mean 
value of the friction coefficient ( in the range (0.44 < e < 1.0) increases from about 0.15 to 
1.35 . Figure 3 shows the gradual expansion of the hysteretic limit cycle as the maximum 
temperature gradient is increased from 0.44 to 1.00 . 

Figures 4 and 5 show a bit of the complexity associated with smaller values of the tem- 
perature gradient. Using an initial momentum of unity gives more regular attractors, of 
the type shown to the left. On the other hand, much higher initial momenta give chaotic 
distributions like those shown to the right. This complexity is no doubt related to that seen 
without any temperature gradient at all 2 ^. In that latter case the phase-space distribution is 
divided into an infinite number of coexisting distributions, whose union is Gibbs' canonical 
distribution , 

f _ „-[ q 2 +p 2 +C 2 ]/2 

/Gibbs — e ■ 

The Figures show projections of a strange attractor that forms with e = 0.40 . The 
Lyapunov spectrum in this case is nearly symmetric, so that it is difficult to compute an 
accurate information dimension of the attractor. 
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p(left) = 1 ; p(right) = 5 
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8x8 Momentum/Friction Plots 

FIG. 4: for the nonequilibrium Nose-Hoover oscillator (e = 0.4) is plotted between the limits 

±4 for two different initial conditions. For p = 5 the Lyapunov exponents are roughly ±0.0025 and 
0. For p = 1 the exponents are much larger in magnitude, ±0.0085 and 0. 



p(left) = 1 ; p(right) = 5 




8x8 (q,p) Phase-Plane Plots 

FIG. 5: p(q) is plotted between the limits ±4 for the two different initial conditions of Figure 4 . 
Note the preference of the oscillator for the lower-temperature states to the left of the origin. 

Fortunately, the complex dynamics of the thermostated oscillator can be greatly simpli- 
fied by adding another control variable, a friction coefficient controlling the fourth velocity 
moment^ : 

{ Q = P ; V = -Q - (P - & 3 ; C = V 2 - T(q) ; £ = p 4 - 3p 2 T(q) ; T(q) = 1 + etanh(g) } . 

At equilibrium the extra control variable allows the oscillator to sample the complete canon- 
ical distribution. This works at nonequilibrium too. Figure 6 compares the distributions of 
the two friction coefficients ((, £) for e equal to 0.5 and 1.0 . Even in the latter case the chaos 
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Temperature Gradients 0.5 and 1.0 




10x10 x(z) friction coefficients 

FIG. 6: Friction coefficient distribution for two values of the maximum temperature gradient, 
0.5 and 1.0 . This doubly-thermostated oscillator covers the complete canonical distribution in the 
equilibrium case. Here there is no dependence of the attractor on the initial value of the momentum. 
100,000 points are printed taken from the last half of a 200 million timestep run. The timestep At 
is 0.0002 . 

induced by the two coefficients is enough to prevent collapse of the dynamics onto a limit 
cycle. Although counterintuitive, it appears to be true that a four- dimensional attractor is 
actually much simpler than its three-dimensional counterpart. 



A. Larger Systems and Thermodynamics 

Larger systems fit the pattern to which the small systems hint. The phase-space dis- 
tribution shrinks to a strange attractor. In a system with several thermostated degrees of 
freedom Liouville's Theorem gives the details of the shrinkage^^S : 

(dlnf/dt) = -(<g>/<g>) = J2C = exp[ (S/k) ] . 

Here S is the external entropy production, the heat extracted from the controlled system by 
the thermostats, divided by the thermostat temperature. <8> is a small comoving phase vol- 
ume. £g> has three possible evolutions: it can expand; it can shrink; or it can remain the same. 
The last possibility is the equilibrium one, with no net heat transfer to the outside world. 
The first possibility (expansion) is ruled out for steady states, as a continually expanding 
phase volume implies catastrophic instability. Only the possibility of continual shrinkage, 
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dissipation, is left. The accessible phase-space states for a nonequilibrium steady state con- 
tinually decrease in number as the volume shrinks (exponentially fast) toward zero. The 
deterministic time-reversible thermostats make possible a simple geometric interpretation 
of the Second Law of Thermodynamics. Nonequilibrium steady states necessarily collapse 
to a zero-volume strange attractor. Thus nonequilibrium states are vanishingly rare. Any 
attempt to reverse the (time-reversible) dynamics would lead to divergence, with a positive 
Lyapunov sum, and a violation of the Second Law. What happens in fact is that, when 
reversed, the dynamics soon breaks its time symmetry and seeks out again the attractor. 
Time-reversible thermostats have deepened our understanding of the Second Law-i>2Z. 

VI. SUMMARY 

The paradoxical reversibility properties of Newtonian and Hamiltonian mechanics can 
be modeled with bit-reversible algorithms. Such algorithms don't exist in cases where the 
phase volume changes, where the mechanics is thermostated. In the latter case Lyapunov 
instability seeks out the unstable strange attractor, more stable still than is its repeller twin, 
leading to a simple geometric understanding of the Second Law of Thermodynamics for open 
systems. 

The symmetry breaking revealed by strong Shockwaves suggests that a deepened under- 
standing of isolated systems can come from study of the local Lyapunov spectrum. Both of 
these problem areas, nonequilibrium conservative systems and nonequilibrium open systems, 
suggest many interesting research opportunities for the future. 
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